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Abstract 


The  resampling  of  a  signal  involves  the  conversion  from  the  initial  sampling  rate  to  a 
new  and  different  one.  This  is  often  necessary  in  practical  applications  because  the 
sampling  rate  is  often  fixed  while  the  desired  sampling  rate  may  depend  on  the 
application  or  a  signal  parameter,  such  as  the  symbol  rate  or  channel  spacing.  This 
problem  often  arises  in  software  radios  systems  where  the  sampling  rate  is  determined 
by  hardware  constraints.  Although  resampling  algorithms  have  been  developed  and 
implemented  in  commercial  software  libraries,  such  as  the  Intel  Signal  Processing 
Library  [2],  these  algorithms  often  have  undesirable  limitations.  For  example,  a  typical 
constraint  is  that  the  data  blocks  processed  must  be  an  integer  multiple  of  the 
downsampling  rate.  This  restriction  simplifies  the  indexing  in  the  code  and  reduces  the 
complexity  of  the  implementation.  However,  it  decreases  the  flexibility  and  usefulness 
of  the  resulting  program,  since  in  some  applications,  the  length  of  the  input  data  blocks 
may  not  be  an  integer  multiple  of  the  downsampling  rate.  For  such  applications,  there  is 
a  need  for  a  resampling  program  that  imposes  no  restriction  on  the  length  of  the  input 
data  blocks.  The  MATLAB  built-in  function  UPFIRDN.dll  (or  RESAMPLE.m)  is 
designed  for  resampling  only  one  block  of  input  data.  It  does  not  provide  information 
on  the  state  variables  of  the  filters  and,  therefore,  cannot  be  used  to  resample  a  very 
large  data  file. 

This  report  describes  a  fairly  general  resampling  algorithm  useful  in  many  practical 
signal  processing  applications.  To  obtain  the  desired  flexibility,  we  have  implemented 
the  periodically  time  varying  FIR  filter  structure  in  such  a  way  that  it  keeps  track  of  the 
state  variables  and  also  makes  provisions  for  input  blocks  of  arbitrary  length.  In 
addition  to  its  flexibility,  the  algorithm  achieves  a  relatively  high  processing  throughput 
when  implemented  in  software. 
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Resume 


Le  reechantillonnage  d’un  signal  implique  la  conversion  du  taux  d’echantillonnage 
initial  a  un  taux  different.  C’est  souvent  une  necessite  dans  certaines  applications 
pratiques  dont  le  taux  d’echantillonnage  est  fixe,  alors  qu’il  faudrait  que  ce  taux 
d’echantillonnage  varie  selon  l’application  ou  un  parametre  du  signal,  tels  le  debit  des 
symboles  ou  1’espacement  entre.  C’est  un  probleme  qui  surgit  ffequemment  pour  les 
plates-formes  radio  logicielle  pour  lesquelles  le  taux  d’echantillonnage  est  determine 
par  des  contraintes  materielles.  Bien  qu’on  ait  developpe  des  algorithmes  de 
reechantillonnage  dans  le  commerce  et  qu’on  les  ait  mis  en  application  dans  des 
bibliotheques  de  logiciels,  comme  la  Intel  Signal  Processing  Library,  ces  algorithmes 
sont  souvent  a  la  merci  de  limitations  indesirables.  Une  contrainte  frequente  est 
1’obligation  que  les  blocs  de  donnees  soient  un  multiple  entier  du  taux  de 
sous-echantillonnage.  Une  telle  restriction  facilite  1’ indexation  du  code  et  reduit  la 
complexite  de  la  mise  en  application.  Cependant,  le  programme  resultant  est  alors 
moins  polyvalent  et  moins  utile,  car  dans  certaines  applications,  la  longueur  des  blocs 
de  donnees  d’ entree  ne  peut  pas  etre  un  multiple  entier  du  taux  de 
sous-echantillonnage.  Pour  ces  applications,  il  faut  recourir  un  programme  de 
reechantillonnage  qui  n’ impose  aucune  restriction  quant  a  la  longueur  des  blocs  de 
donnees  d’ entree.  La  fonction  integree  de  MATLAB  UPFIRDN.dll  (ou 
RESAMPLE.m)  a  ete  ecrite  pour  ne  reechantillonner  qu’un  seul  bloc  d’entree  de 
donnee.  Elle  ne  fournit  aucune  information  sur  les  variables  d’etat  des  filtres  et,  par 
consequent,  elle  ne  peut  servir  a  reechantillonner  un  tres  grand  fichier  de  donnees. 

Le  present  rapport  decrit  un  algorithme  de  reechantillonnage  suffisamment  general  qui 
pourra  etre  utile  dans  de  nombreuses  applications  pratiques  de  traitement  des  signaux. 
Pour  obtenir  la  souplesse  souhaitee,  nous  avons  applique  la  structure  de  filtre  FIR 
variable  de  telle  maniere  qu’elle  controle  les  variables  d’etat  et  permet  des  blocs 
d’entree  de  longueur  quelconque.  En  plus  de  cette  souplesse,  1’ application  logicielle  de 
cet  algorithme  atteint  un  debit  de  traitement  relativement  eleve. 


Executive  summary 


The  resampling  of  a  signal  involves  the  conversion  from  the  initial  sampling  rate  to  a 
new  and  different  one.  This  is  necessary  in  many  practical  applications  because  the 
sampling  rate  is  often  fixed  while  the  desired  sampling  rate  may  depend  on  the 
application  or  a  signal  parameter,  such  as  the  symbol  rate  or  channel  spacing.  It  is  a 
particular  issue  in  digital  filter  bank  receivers  used  to  detect  the  presence  of  narrowband 
signals  occupying  channels  over  a  range  of  frequencies.  To  minimize  computational 
cost,  filter  bank  receivers  often  use  frequency  domain  techniques  based  on  the  Fast 
Fourier  Transform  (FFT)  algorithm.  Since  it  is  usually  desirable  to  select  FFTs  to  have 
a  length  that  is  a  power  of  two,  the  channelization  frequencies  may  not  align  with  the 
channel  spacings  used  in  a  certain  frequency  bands.  The  situation  is  complicated 
further  by  the  introduction  of  new  channel  spacing  standards,  e.g.,  8.33  kHz.  An 
example  of  this  problem  occurs  with  the  Agilent  Blackbird  signal  analysis  system.  For 
the  standard  sampling  rate  of  20.48  megasamples/s,  FFT  bin  sizes  are  constrained  to  the 
set ...  1.25  kHz,  2.5  kHz,  5  kHz,  10  kHz, ...  with  the  result  that  a  channelization  such  as 
8.33  kHz  cannot  be  accommodated  directly. 

Although  resampling  algorithms  have  been  developed  and  implemented  in  commercial 
software  libraries,  such  as  the  Intel  Signal  Processing  Library  [2],  these  algorithms 
often  have  undesirable  limitations.  For  example,  a  typical  constraint  is  that  the  data 
blocks  processed  must  be  an  integer  multiple  of  the  downsampling  rate.  This  restriction 
simplifies  the  indexing  in  the  code  and  reduces  the  complexity  of  the  implementation. 
However,  it  decreases  the  flexibility  and  usefulness  of  the  resulting  program,  since  in 
some  applications,  the  length  of  the  input  data  blocks  may  not  be  an  integer  multiple  of 
the  downsampling  rate.  Another  implementation  of  a  resampling  algorithm,  the 
MATLAB  built-in  function  UPFIRDN.dll  (or  RESAMPLE.m)  is  designed  for 
resampling  only  single  blocks  of  input  data.  It  does  not  provide  information  on  the  state 
variables  of  the  filters  and,  therefore,  cannot  be  used  to  resample  a  very  large  data  file. 
Consequently,  there  is  a  need  for  a  resampling  program  that  imposes  no  restriction  on 
the  length  of  the  input  data  blocks.  To  obtain  a  more  flexible  resampling  program  we 
have  implemented  the  periodically  time  varying  FIR  filter  structure  in  such  a  way  that  it 
keeps  track  of  the  state  variables  and  also  makes  provisions  for  input  blocks  of  arbitrary 
length.  The  implementation  is  carried  out  in  MEX  in  the  MATLAB  environment  with 
the  core  of  the  program  written  in  C.  In  addition  to  utilizing  the  very  efficient 
periodically  time- varying  FIR  filter  structure,  an  attempt  is  made  to  minimize  data 
movement  in  the  implementation.  It  turns  out  that  the  implementation  is  not  only  very 
flexible  but  also  works  approximately  twice  as  fast  as  the  Matlab  built-in  function 
UPFIRDN.dll  for  non-trivial  resampling  ratios  (that  is,  L  >  1  and  M  >  1,  with  L 
denoting  the  upsampling  factor  and  M  denoting  the  downsampling  factor). 

The  new  algorithm  is  particularly  useful  in  digital  receiver  systems  where  a  continuous 
stream  of  signal  data  must  be  processed,  for  example,  in  the  recovery  of  message 
content  from  a  signal  of  arbitrary  duration. 
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Sommaire 


Le  reechantillonnage  d’un  signal  implique  la  conversion  du  taux  d’echantillonnage 
initial  a  un  taux  different.  C’est  souvent  une  necessite  dans  les  applications  pratiques 
dont  le  taux  d’echantillonnage  est  fixe,  alors  qu’il  faudrait  que  ce  taux 
d’echantillonnage  varie  selon  l’application  ou  un  paramfctre  du  signal,  tel  le  debit  des 
symboles  ou  l’espacement  entre.  C’est  un  probleme  bien  precis  dans  les  recepteurs  a 
banque  de  filtres  numeriques  que  Ton  utilise  pour  detecter  la  presence  de  signaux  a 
largeur  de  bande  etroite  qui  occupent  des  dans  une  plage  de  frequences.  Pour  reduire  au 
minimum  le  cout  des  calculs.  les  recepteurs  a  banque  de  filtres  font  souvent  appel  des 
techniques  du  domaine  de  la  frequence  fondees  sur  l’algorithme  de  la  transformee 
rapide  de  Fourier  (TRF).  Comme  il  est  generalement  souhaitable  de  selectionner  des 
TRF  pour  avoir  une  longueur  qui  soit  une  puissance  de  2,  la  resoloution  en  frequence 
de  decoupage  peut  ne  pas  correspondre  avec  l’espacement  tel  des  voies  de  certaines 
bandes  de  frequences.  La  situation  se  complique  du  fait  de  1’ introduction  de  nouvelles 
normes  d’ecartement  de  voies,  tel  8,33  kHz.  Un  exemple  de  ce  probleme  se  retrouve 
dans  le  systeme  d’ analyse  de  signaux  de  Agilent  Blackbird.  Pour  les  taux 
d’echantillonnage  standard  de  20.48  megaechantillons/s,  les  longueurs  des  intervalles 
de  la  TRF  sont  restreintes  a  1.25  kHz,  2.5  kHz,  5  kHz,  10  kHz,  de  sorte  qu’un 
decoupage  de  voies  comme  8.33  kHz  ne  peut  etre  decrit  directement. 

Bien  qu’on  ait  developp6  des  algorithmes  de  reechantillonnage  dans  le  commerce  et 
qu’on  les  ait  mis  en  application  dans  des  bibliotheques  de  logiciels,  comme  la  Intel 
Signal  Processing  Library,  ces  algorithmes  sont  souvent  a  la  merci  de  limitations 
indesirables.  Une  contrainte  frequente  est  l’obligation  que  les  blocs  de  donnees  soit  un 
multiple  entier  du  taux  de  sous-echantillonnage.  Une  telle  restriction  facilite 
1’ indexation  du  code  et  reduit  la  complexity  de  la  mise  en  application.  Cependant,  le 
programme  resultant  est  alors  moins  polyvalent  et  moins  utile,  car  dans  certaines 
applications,  la  longueur  des  blocs  de  donnees  d’entree  ne  peut  pas  etre  un  multiple 
entier  du  taux  de  sous-echantillonnage.  La  fonction  integree  de  MATLAB 
UPFIRDN.dll  (ou  RESAMPLE.m)  a  ete  ecrite  pour  ne  reechantillonner  qu’un  seul  bloc 
d’entree  de  donnee.  Elle  ne  fournit  aucune  information  sur  les  variables  d’etat  des 
filtres  et,  par  consequent,  elle  ne  peut  servir  reechantillonner  un  tres  grand  fichier  de 
donnees.  Un  programme  de  reechantillonnage  qui  n’impose  aucune  restriction  quant  la 
longueur  des  blocs  d’entree  est  done  justifie.  Pour  realiser  tel  programme  plus 
polyvalent,  nous  avons  applique  une  structure  de  filtre  HR  temporelle  a  variation 
periodique  qui  controle  les  variables  d’etat  et  permet  des  blocs  d’entree  de  longueur 
variable.  La  mise  en  application  est  realisee  en  MEX  dans  un  environnement 
MATLAB,  le  noyau  du  programme  etant  redige  en  C.  En  plus  d’exploiter  la  structure 
de  filtre  FIR  temporelle  a  variation  periodique  tres  performante,  nous  avons  tente  de 
reduire  au  minimum  le  mouvement  des  donnees  dans  cette  mise  en  application.  Cette 
mise  en  application  s’ est  revelee  moins  souple  que  prevu,  mais  quand  meme  deux  fois 
plus  rapide  que  la  fonction  integree  MatLab  UPFIRDN.dll  pour  des  rapports  de 
reechantillonnage  non  triviaux  (soit  L  >  1  et  M  >  1  o  L  est  le  facteur  de 
sur-echantillonnage  et  M,  le  facteur  de  sous-echantillonnage). 
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Ce  nouvel  algorithme  est  particulierement  utile  pour  les  recepteurs  numeriques  dans 
lesquels  le  flot  continu  de  donnees  numerique  doit  etre  traite,  par  exemple,  lors  de  la 
recuperation  du  contenu  d’un  message  dans  un  signal  de  duree  quelconque. 
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1.  INTRODUCTION 


The  resampling  of  a  signal  involves  the  conversion  from  the  initial  sampling  rate  to  a 
new  and  different  one.  This  is  often  necessary  in  practical  applications  because  the 
sampling  rate  is  often  fixed  while  the  desired  sampling  rate  may  depend  on  the 
application  or  a  signal  parameter,  such  as  the  symbol  rate  or  channel  spacing.  It  is  a 
particular  issue  in  digital  filter  bank  receivers  used  to  detect  the  presence  of  narrowband 
signals  occupying  channels  over  a  range  of  frequencies.  To  minimize  computational 
cost,  filter  bank  receivers  often  use  frequency  domain  techniques  based  on  the  Fast 
Fourier  Transform  (FFT)  algorithm.  Since  it  is  usually  desirable  to  select  FFTs  to  have 
a  length  that  is  a  power  of  two,  the  channelization  frequencies  may  not  align  with  the 
channel  spacings  used  in  a  certain  frequency  bands.  The  situation  is  complicated 
further  by  the  introduction  of  new  channel  spacing  standards,  e.g.,  8.33  kHz.  An 
example  of  this  problem  occurs  with  the  Agilent  Blackbird  signal  analysis  system.  For 
the  standard  sampling  rate  of  20.48  megasamples/s,  FFT  bin  sizes  are  constrained  to  the 
set ...  1.25  kHz,  2.5  kHz,  5  kHz,  10  kHz, ...  with  the  result  that  a  channelization  such  as 
8.33  kHz  cannot  be  accommodated  directly. 

Although  resampling  algorithms  have  been  developed  and  implemented  in  commercial 
software  libraries,  such  as  the  Intel  Signal  Processing  Library  [2],  these  algorithms 
often  have  undesirable  limitations.  For  example,  a  typical  constraint  is  that  the  data 
blocks  processed  must  be  an  integer  multiple  of  the  downsampling  rate.  This  restriction 
simplifies  the  indexing  in  the  code  and  reduces  the  complexity  of  the  implementation. 
However,  it  decreases  the  flexibility  and  usefulness  of  the  resulting  program,  since  in 
some  applications,  the  length  of  the  input  data  blocks  may  not  be  an  integer  multiple  of 
the  downsampling  rate.  Another  implementation  of  a  resampling  algorithm,  the 
MATLAB  built-in  function  UPFIRDN.dll  (or  RESAMPLE.m)  is  designed  for 
resampling  only  single  blocks  of  input  data.  It  does  not  provide  information  on  the  state 
variables  of  the  filters  and,  therefore,  cannot  be  used  to  resample  a  very  large  data  file. 
Consequently,  there  is  a  need  for  a  resampling  program  that  imposes  no  restriction  on 
the  length  of  the  input  data  blocks. 

To  obtain  a  more  flexible  resampling  program  than  the  ones  currently  available,  we 
have  implemented  the  periodically  time  varying  FIR  filter  structure  in  such  a  way  that  it 
keeps  track  of  the  state  variables  and  also  makes  provisions  for  input  blocks  of  arbitrary 
length.  The  implementation  is  carried  out  in  MEX  in  the  MATLAB  environment  with 
the  core  of  the  program  written  in  C.  In  addition  to  utilizing  the  very  efficient 
periodically  time- varying  FIR  filter  structure,  an  attempt  is  made  to  minimize  data 
movement  in  the  implementation.  It  turns  out  that  the  implementation  is  not  only  very 
flexible  but  also  works  approximately  twice  as  fast  as  the  Matlab  built-in  function 
UPFIRDN.dll  for  non-trivial  resampling  ratios  (that  is,  L  >  1  and  M  >  1,  with  L 
denoting  the  upsampling  factor  and  M  denoting  the  downsampling  factor). 

This  report  is  organized  as  follows.  Section  2  summarizes  the  theoretical  basis  of 
digital  sampling  rate  conversion  and  derives  the  periodically  time-varying  FIR  filter 
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structure.  Section  3  discusses  the  special  implementation  of  the  periodically 
time- varying  FIR  filter  structure  in  MEX  (named  UPDNC6.c).  Section  4  demonstrates 
the  usage  of  two  MATLAB  functions  UPDNM6.m  and  RESAM6.m  based  on 
UPDNC6.dll.  Section  5  presents  a  comparison  of  the  performance  of  our 
implementation  UPDNC6.dll  with  that  of  the  MATLAB  implementation  UPFIRDn.dll. 
Finally,  section  6  presents  some  concluding  remarks. 

2.  BASIC  CONCEPTS  OF  DIGITAL  SAMPLING  RATE 
CONVERSION 


Let  L  >  0  and  M  >  0  be  two  arbitrary  positive  integers  which  are  relatively  prime,  that 
is,  the  greatest  common  divisor  of  L  and  M  is  1.  A  baseband  digital  signal  x(n),  with 
sampling  rate  Fo ,  can  be  upsampled  by  a  factor  of  L  and  then  downsampled  by  a  factor 
of  M  to  obtain  another  digital  signal  y(m)  with  sampling  rate  Fi  =  F0  (^).  In  theory, 
upsampling  (also  called  interpolation)  and  downsampling  (also  called  decimation)  are 
two  separate  processes  and  upsampling  must  precede  downsampling  in  order  to 
preserve  the  spectrum  of  the  signal.  In  practice,  the  upsampling  and  downsampling 
processes  are  combined  as  one  and  can  be  implemented  efficiently  via  the  periodically 
time- varying  FIR  filter  structure  as  discussed  in  the  classical  survey  paper  [1].  In  this 
section,  the  main  concepts  of  digital  sampling  rate  conversion  are  discussed  and  the 
periodically  time- varying  FIR  filter  structure  is  derived. 


2.1  Downsampling  by  an  Integer  Factor  of  M. 


Let  x(n)  be  a  baseband  digital  signal  with  sampling  rate  Fo.  The  sampling  rate  of  x(n) 
can  be  decreased  by  a  factor  of  M  by  retaining  only  one  sample  in  every  M  samples  in 
x(n).  Specifically,  if  the  sequence  y(m)  is  defined  by  y(m)  =  x(mM),  we  say  that 
y{m)  is  obtained  by  downsampling  x(n)  by  a  factor  of  M  (or  by  decimating  x(n)  by  a 
factor  M).  Direct  decimation  of  the  signal  x(n)  without  first  passing  it  through  an 
appropriately  designed  low-pass  filter  will  in  general  lead  to  aliasing  in  the  resultant 
signal  y(m).  To  demonstrate  this,  one  can  calculate  the  z  transform  of  y(m).  Let  the 
sequence  x'(n)  be  defined  by 


(1) 


x'(n) 


f  x(n),  n  is  an  integer  multiple  of  M, 
0,  otherwise 


It  can  be  verified  that 


(2) 


M— 1 


An)  =  x(n)  \  ±  Y,  ^ ^ 


1=0 


and 


(3) 


y(m)  =  x  (mM) 
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Hence,  Y ( z ),  the  z  transform  of  y(m),  can  be  calculated  as  follows: 

oo 

y(z)  =  Y  y(m)^_m 


(4) 


=  Y  x'{mM)z-m 

771=— OO 
OO 

-  Y  x'(m)z~m/M 


771=  — OO 
OO 


M-l 


\ 


-  E+)  iE^" 

771=  —  OO  L  i  =  0  J 


-m/M 


M-l 


1=0 

M-l 


-  iW  —  i 

4E  E  x{Tn)e?2vlm/M  z~m/M 


L771=  —  OO 


1  iva  —  ± 

=  X{e~j2nl/Mz1/M) 


1=0 


where  X{z)  is  the  z  transform  of  x(ri).  Evaluating  Y (z)  on  the  unit  circle,  z  =  eju, 
uj  e  [-7T,  7 r],  yields  the  Fourier  transform  of  y(m),  i.e., 

M-l 

(5)  Y(ej“)  =  j;Y  X(e^~2^/M) 

1=0 


^Frorn  Eq.  5  it  can  be  seen  that  in  general  y(m)  is  an  aliased  version  of  the  original 
signal  x(n). 


To  avoid  aliasing  in  y(m),  it  is  necessary  to  first  low-pass  filter  x(n)  and  then  perform 
decimation.  In  fact,  let  hd{n)  be  the  impulse  response  of  a  low-pass  FIR  filter(called 
the  prototype  downsampling  filter)  and  Hd(z)  be  its  z  transform.  Let  the  signal  w(n)  be 
the  output  of  the  filter  Hd(z)  with  input  x(n)  and  let  y(m)  be  the  signal  obtained  by 
decimating  w(n)  by  a  factor  of  M  (that  is,  y{m)  —  w(Mm)).  Let  X(z),  W (z)  and 
Y(z)  be  the  z  transforms  of  x(n),  w(n)  and  y(rrt)  respectively.  It  follows  that 

(6)  W(z)  -  Hd(z)X(z) 

and  (see  Eq.  4,  with  W(z )  replacing  X(z )) 

M-l 

Y(z)  =  jjYl  W{e~jM/Mzl/M) 

1=0 

M-l 

Y(z)  =  £  W^e~jM/Mzl/M) 

1=0 
-  M-l 

=  —Y  Hd(e-^Mzl/M)X(e-^Mz^M) 

^  1=0 


(7) 

Hence 

(8) 
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Figure  1:  Downsampler  by  a  factor  of  M  and  prototype  lowpass  filter. 


Evaluating  Y ( z )  on  the  unit  circle  z  =  e^,  u>  e  [— w,  7r],  yields  the  result 
1  M—l 

Y(+)  =  jjYl  Hd(ej{“~2T,l)/M)X(ej{bJ-27rl)/M) 
l=o 

(9)  =  T  ^Hd(eju/M)X{^M)  +  Hd(e*<a-**VM )X  (e^"2*^)  +  ■■■ 

If  the  prototype  downsampling  filter  Hd(z)  is  designed  in  such  a  way  that 


(10) 


Hits* 


W\  J  1 1 
j-\  0, 


Ml  < 

X?  <  Ml  <  7T 


then  the  terms  with  l  ^  0  in  Eq.  9  are  removed  and  Y (e-7^)  becomes 

(ii) 

Hence  by  passing  x(n)  through  an  appropriately  designed  FIR  prototype  downsampling 
filter  Hd(z)  and  then  performing  decimation,  aliasing  in  y(m)  is  eliminated.  The 
process  of  downsampling  a  digital  signal  by  a  factor  of  M  is  illustrated  in  Fig.  L 


The  time  domain  relationship  among  x(n),  y(m),  hd(k)  is 


oo 

(12)  y(m)  =  w(mM)  =  ^  hd(k)x(mM  —  k) 

k=  —  oo 
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2.2  Upsampling  by  an  Integer  Factor  of  L. 


Upsampling  is  the  dual  operation  of  downsampling.  Suppose  it  is  desired  to  increase 
the  sampling  rate  of  x(n)  by  an  integer  factor  of  L.  This  can  be  accomplished  by 
inserting  L  -  1  zeros  between  each  pair  of  samples  in  x(n)  and  then  low-pass  filtering 
the  resultant  sequence.  In  fact,  let  the  digital  signal  w(n)  be  defined  by 


(13) 


x(^)}  n  is  an  integer  multiple  of  L, 
0,  otherwise 


The  signal  w(ri)  is  obtained  by  inserting  L  -  1  zeros  between  each  pair  of  samples  of 
x(n).  Its  z  transform  W(z)  is  given  by 

oo 

W(z )  =  ^  w(n)z~n 


OO 

=  x(n)z~nL 
71—  —  OC 

(14)  =  X(zL ) 

where  X  ( z )  is  the  z  transform  of  x(n).  Evaluating  W (z)  on  the  unit  circle  ,  z  =  eJW, 
uj  €  [ — 7r,  7r],  gives  the  Fourier  transform  of  w(n) 

(15)  W{e>“)  -  X(ejuL) 

Clearly,  the  spectrum  of  w(n )  preserves  that  of  the  original  signal  x(n)  on  the  interval 
[0,  j]  but  also  contains  harmonic  images  of  the  spectrum  of  x(n)  on  the  intervals 
[f  ’  Z 7rl’  ‘ >  [^r^7r.  7T}<  which  must  be  removed.  This  is  illustrated  in  Figure  2  for  the 
case  L  =  4. 


To  eliminate  the  unwanted  harmonic  images  of  the  spectrum  of  x(n),  it  is  necessary  to 
filter  the  signal  w(n)  with  a  low-pass  filter  Hu(z)  (called  the  prototype  upsampling 
filter)  which  approximates  the  ideal  frequency  response 


(16) 


Hu(ejn 


M  <  f’ 

o,  f  <  M  <7T 


where  G  =  L  is  a  necessary  scaling  factor  (see  Eq.  33  [1]).  The  process  of  upsampling 
a  digital  signal  x(n)  by  a  factor  of  L  is  illustrated  in  Fig.  3. 


The  time  domain  relationship  among  x(n),  y(m),  hu(k),  where  hu(k)  is  the  impulse 
response  of  the  prototype  upsampling  filter  Hu(z),  is  given  by 

OO 

y{m )  =  ^  hu(m  -  k)w(k) 

k=- oo 

=  ^  hu(m  -  k)x(k/L) 

i  is  an  integer 

OO 

(17)  =  hu(m  -  kL)x(k) 

k=— oo 
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Figure  3:  Upsampler  by  factor  of  L  followed  by  lowpass  filter  to  remove  harmonics. 
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2.3  Resampling  by  a  Rational  Factor  of  jj. 

Changing  the  sampling  rate  of  a  digital  signal  by  an  arbitrary  rational  factor  of  , 
where  L  and  M  are  relatively  prime  positive  integers,  is  accomplished  by  first 
increasing  the  sampling  rate  by  a  factor  of  L  and  then  decreasing  the  sampling  rate  of 
the  resultant  signal  by  a  factor  of  M.  As  demonstrated  in  the  preceding  two 
subsections,  for  the  upsampling  stage,  the  prototype  upsampling  filter  Hu{z)  should 
approximate  the  frequency  response  characteristic 


(18) 


M  <  b 


and  for  the  downsampling  stage,  the  prototype  downsampling  filter  Hd(z)  should 
approximate  the  frequency  response  characteristic 


(19) 


M  -  M  ’ 

If  <  M  <  7T 


The  complete  resampling  process  is  carried  out  by  first  inserting  L  -  1  zeros  between 
each  pair  of  samples  of  the  input  signal  x(n),  low-pass  filtering  the  resultant  signal  by 
the  filter  Hu(z),  then  low-pass  filtering  the  output  signal  of  Hu(z)  by  the  filter  Hd(z) 
and  finally  decimating  the  output  signal  of  Hd(z)  by  a  factor  of  M.  This  is 
demonstrated  in  Fig.  4.  This  process  can  be  further  simplified  by  replacing  the  cascade 
of  the  prototype  filters  Hu(z )  and  Hd(z)  by  one  low-pass  filter  H(z)  —  Hu(z)Hd(z), 
which  has  the  frequency  response  characteristic 


(20) 


H(z)  =  Hu(z)Hd(z)  S* 


L'  M  ^  ma x(M,L)  > 

°>  ma x(M,i)  <  lwl  -  7r 


The  low-pass  filter  H(z)  is  called  the  prototype  resampling  filter.  Conceptually, 
resampling  is  now  reduced  to  three  steps:  first  padding  L  —  1  zeros  between  each  pair 
of  samples  of  the  input  signal  x(n),  then  low-pass  filtering  the  resultant  signal  by  the 
resampling  filter  H(z),  and  finally  decimating  the  output  signal  of  H(z)  by  a  factor  M. 
This  is  demonstrated  in  Fig.  5. 


Figure  4:  The  process  of  resampling  a  signal  by  a  factor  of  L/M  by  upsampling,  filtering,  and  downsampling. 

Let  h(k)  be  the  impulse  response  of  the  prototype  resampling  filter  H{z)  defined  by 
Eq.  20.  The  time  domain  relationship  among  x(n),  h(k)  and  y(m)  is 

(21)  y{m)  =  v{mM) 
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x(n) 

H 

w(m) 

Low  Pass 

Filter 

H(z) 

Decimator 

y(m)  “  v(mM) 

Interpolator 

( 

) 

s 

_ i 

L _ 

-7i  -rc/max(L,M)  0  jc/max(L.M)  n 


Figure  5:  The  process  of  resampling  a  signal  with  a  single  prototype  lowpass  filter. 
where  v(n)  is  the  output  of  the  filter  H(z)  and 

oo 

(22)  v(m)  =  ^  h(m  —  kL)x(k) 

fc=— oo 

Therefore 

oo 

(23)  y(m)=  ]T  h(mM-kL)x(k) 

k=- oo 

Equation  23  is  the  mathematical  basis  for  implementing  digital  sampling  rate 
conversion. 


2.4  The  Periodically  Time-varying  FIR  Filter  Structure  for 
Implementing  Digital  Sampling  Rate  Conversion. 


Equation  23  can  be  put  in  a  form  more  amenable  to  hardware  or  software 
implementation.  Making  the  change  of  variables 

.  ( mM\ 

(24)  k  =  floor  (  —j—  j  -  n 

where  floor  (m^-)  denotes  the  largest  integer  less  than  or  equal  to  the  rational  number 
one  can  write 


(25) 


mM  —  kL  —  mM  -  ^floor  —  nj  L 

=  nL  +  mM-  ^floor  L 


—  nL  4-  mM  ©  L 
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where  mM  ©  L  =  mM  -  (floor  (2^))  L  denotes  the  remainder  of  mM  after  being 
divided  by  L  (value  of  mM  modulo  L).  Substituting  Eqs.  24  and  25  in  23,  we  obtain 

(26)  y{m)  =  ^  h(nL  +  mM  ©  L)x  f  floor  \  ~ 

n=-oo  \  \  /  / 


This  is  the  actual  resampling  equation  that  is  implemented  in  hardware  or  software. 

In  the  practical  implementation  of  Eq.  26,  the  prototype  resampling  filter  h(k)  is 
assumed  to  be  a  linear  phase  FIR  filter  with  length  LQ  ,  where  Q  >  1  is  a  positive 
integer.  That  is,  the  length  of  h(k)  is  constrained  to  be  an  integer  multiple  of  the 
upsampling  rate  L.  Under  this  constraint,  the  algorithm  Eq.  26  leads  naturally  to  the 
periodically  time-varying  FIR  filter  structure  for  sampling  rate  conversion. 

Some  new  notation  will  now  be  introduced.  The  prototype  resampling  filter  h(k), 
where  0  <  k  <  LQ  —  1,  can  be  partitioned  into  L  polyphase  filters  hm(n) 

(0  <  m  <  L  -  1),  each  of  length  Q,  i.e., 


(27) 


hm(n)  =  h(nL  -I-  m),  0  <  n  <  Q  —  1 


In  other  words,  hm(n)  is  obtained  by  decimating  the  sequence  h(k)  by  a  factor  of  L. 
To  visualize  the  relationship  between  the  prototype  filter  h(k)  (with  length  N  =  LQ) 
and  the  polyphase  filters  hm[n)  (each  with  length  Q),  the  Q  x  L  matrix  H  is 
constructed,  i.e.. 


(  HO) 

Ml) 

M2) 

•  •  •  h(L  -  1)  \ 

h(L) 

h(L  +  1) 

h(L  +  2) 

•••  h(2L  —  1) 

H  = 

h(nL) 

h(nL  1) 

h(nL  +  2) 

•••  h((n  +  1)L  —  1) 

h((Q  -  1  )L) 

h((Q-l)L  +  l) 

h((Q  -  1)L  +  2) 

•••  h(QL  -  1)  / 

(28) 


H  has  L  columns,  the  first  column  consisting  of  the  taps  of  the  polyphase  filter  ho(n), 
the  second  column  consisting  of  the  taps  of  the  polyphase  filter  hi(n)  and  so  on.  On  the 
other  hand,  the  filter  taps  h(k)  can  be  recovered  from  H  by  concatenating  the  row 
vectors  of  H  from  the  first  to  the  last  row.  Define  two  integer  sequences  p(m )  and  q(m) 
by  setting 

(29)  p(m)  —  ©  L,  q(m)  =  floor 

It  can  be  verified  that  p(m)  is  a  periodic  sequence  with  period  L  and  q(m)  has  the 
simple  property  that  q{m  +  kL)  =  kM  +  q[m),  i.e., 

(30)  p{m  +  kL)=p(m)=p(m®L),  q(m  +  kL)  =  kM  +  q(m), 
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Using  this  new  notation,  Eq.  26  can  be  rewritten  as 

y(m)  =  ^  h(nL  +  mM  0  L)x  (  floor  [ 

n= 0  k  V  L 

Q- 1 

=  hp{m)(n)x(q(m)  -  n) 

Tl~  0 

Q-i 

(31)  =  hp{mS>L)(n)x(q(m)  -  n) 

n— 0 

Let 

(32)  gm{n)  =  hp(mQL)(n),  1  <m<L,  0  <n  <  Q  -  1 
It  follows  from  Eqs.  31  and  32  that 

<2-1 

(33)  y(kL  +  m)  —  ^  gm(n)x(kM  +  q{m)  -  n) 

n= 0 

where  1  <  m  <  L  and  k  >  0.  The  equation  (33)  can  be  interpreted  as  representing  a 
periodically  time- varying  FIR  filter.  To  illustrate  this,  write  out  Eq.  33  for  1  <  m  <  L, 


k  =  0, 

1,  as  follows: 

y(!) 

= 

pi(0)x(g(l))  +5i(l)x(g(l)  -  1)  + 

- H  9i{Q  -  l)z(g(l)  ~(Q~  1)) 

y(  2) 

= 

g2(0)x{q{2))  +  g2(l)x(q{2)  -  1)  + 

"•+fl2(Q-l)*(«(2)-(Q-l)) 

y(m) 

= 

9m(0)x(q(m ))  +  gm(l)x{q(m)  -  1)  +  •  •  •  +  gm{Q  -  1  )x{q{m)  -  (Q  -  1)) 

y(L) 

= 

gL(0)x(q(L))  +  gL(l)x(q(L)  -!)  +  •••  +  gL(Q  -  1  )x(q(L)  -  (Q  -  1)) 

— 

gL{0)x(M)  +  gL(l)x(M  -  1)  +  •  • 

•  +  9l{Q  -  1)*(M  -  (Q  -  1)) 

y(L  + 1) 

= 

gi(0)x(M  +  q{  1))  +  •  •  •  +  gi{Q  - 

l)x(M  +  q(l)-(Q-l)) 

y{L  +  2) 

= 

g2{0)x(M  +  q( 2))  H - -f  92(Q  — 

l)x(M  +  q{2)  —  (Q  —  1)) 

y{L  +  m) 

r= 

gm(0)x(M  +  q(m))  +  •  •  •  +  gm(Q 

—  l)x(M  +  q(m)  —  (Q  —  1)) 

y(2  L) 

= 

gL(0)x(2M)  +  gL(l)x(2M  -  1)  + 

- h  9l(Q  —  l)x(2M  —  (Q  —  1)) 

(34) 

The  first  equation  in  34  shows  that  y(l)  is  obtained  as  a  weighted  sum  of  Q  sequential 
samples  of  x(n)  starting  at  the  sample  x(q(  1))  and  going  backwards  in  n  sequentially. 
The  weighting  coefficients  are  the  taps  of  the  polyphase  filter  £i(n),  0  <  n  <  Q  -  1. 
Similarly,  y( 2)  is  obtained  as  a  weighted  sum  of  Q  sequential  samples  of  x(n)  starting 
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at  the  sample  x(q( 2))  and  going  backwards  in  n  sequentially.  The  weighting 
coefficients  are  the  taps  of  a  different  polyphase  filter,  namely  g2(n).  This  pattern 
continues  until  the  sample  y(L),  with  y(L)  computed  as  a  weighted  sum  of  Q 
sequential  samples  of  x(n)  starting  at  the  sample  x(q(L))  =  x(M)  and  going 
backwards  in  n  sequentially.  The  weighting  coefficients  in  the  computation  of  y(L)  are 
the  taps  of  the  polyphase  filter  gt(n).  It  is  clear  that  different  sets  of  filter  coefficients 
(namely,  gm(n))  are  used  in  the  computation  of  the  samples  y{m),  m  =  1, 2,  •  •  • ,  L  . 
Now  let  us  look  at  the  computational  procedure  for  the  next  L  samples, 
y(L  +  1),  y(L  +  2),  •  •  • ,  y(2L).  It  can  be  observed  from  the  equations  in  34  that 
y(L  +  1)  is  computed  using  the  same  polyphase  filter  as  y(  1)  (that  is,  fli(n)).  Also 
y(L  +  2)  is  generated  using  the  same  polyphase  filter  as  y(2)  (that  is,  <?2(™))>  etc. 

Hence  the  sample  sequence  y{m)  is  the  output  of  a  periodically  time-varying  FIR  filter 
with  period  L. 

It  is  also  important  to  understand  how  the  input  data  samples  x(n)  enter  into  the 
computation  in  Eq.  33.  Writing  out  Eq.  33  for  y(m)  and  y(m  +  1),  where 
1  <  m  <  L  —  1,  yields 

y(rn)  =  gm(0)x(q(m,))  +  gm(l)x(q(m)  -  1)  + 

(35)  •  •  •  +  gm{Q  ~  l)x{q{m)  -  (Q  -  1)) 

and 

y(m  +  1)  =  g(m+i)(0)x(g(m  +  1))  +  5(m+ i)(l)z(fl(m  +  1)  -  1)  + 

(36)  •  ■  •  +  g(m+i)(Q  -  i)z(g(m  + 1)  -  (Q  - 1)) 

One  can  see  that  y(m)  is  the  weighted  sum  of  Q  sequential  samples  of  x(n)  starting  at 
the  sample  x(q(m))  and  going  backwards  sequentially.  The  samples  from  the  input 
sequence  x(n)  which  are  involved  in  the  weighted  sum  are 

(37)  S  =  [x(g(m)),  x{q(m)  -  1),  •  •  • , x(q(m)  -  (Q  -  2)),  x(q{m)  -  (Q  -  1))] 

The  samples  x(q(m)),  x(q(m)  —  1),  •  ♦  •  ,x(q(m)  —  (Q  —  2)),x(q(m)  —  (Q  —  1))  are 
called  the  state  variables  of  the  time- varying  FIR  filter  in  Eq.  33  and  the  vector  S  is 
called  the  state  variable  buffer  corresponding  to  the  output  sample  y(m).  In  the 
computation  of  y(m  +  1),  the  state  variable  buffer  has  changed  to 

(38)  S  =  [x(q(m  +  1)),  x(q{m  +  1)  -  1),  ■  •  •  ,x(g(m  +  1)  -  (Q  -  1))] 

If  q(m  *F  1)  =  floor  >  q(m)  =  floor  (^),  there  are  q(m  +  1)  -  q(m) 

new  samples  shifted  into  and  q{wn  +  1)  —  <?(m)  old  samples  shifted  out  of  the  state 
variable  buffer  S.  On  the  other  hand,  if  q(m  +  1)  =  q(m),  then  the  state  variable  buffer 
S  remains  unchanged.  It  should  be  noted  that  the  sequence  q(m)  is  non-decreasing  and 
it  is  q(wn)  that  determines  the  input  samples  x(n)  entering  into  the  computation  of  the 
output  sample  y{m).  Recall  that  q(m  +  kL)  =  kM  +  <?(m)  (see  Eq.  30),  hence  for 
l<m<L-lyk>01 

(39)  q{m  +  1  +  kL)  -  q{m  +  kL)  =  kM  +  q(m  +  1)  -  {kM  +  q{m)) 

=  q(m  +  1)  -  q{m) 
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Equation  39  implies  that  the  pattern  in  which  the  data  samples  x(n)  are  shifted  into  and 
out  of  the  state  variable  buffer  S  also  repeats  itself  periodically.  Let  us  examine  in  some 
detail  how  the  state  variable  buffer  S  is  updated  for  the  first  L  samples 
2/(1),  2/(2),  •  •  • ,  y(L).  Denoting  the  contents  of  the  state  variable  buffer  S 
corresponding  to  the  samples  y(l),  y(2),  ■ ,  y(L)  by  Si,  •  •  •,  S^,  respectively,  we  have 

51  =  [x(g(l)),x(9(l)-l),---,a:(g(l)-(Q-l))] 

52  =  [x(q{2)),x{q(2)-l),---,x(q{2)-(Q-l))} 

53  =  Kg(3)),x(g(3)-l),-..,x(g(3)-(Q-l))] 


SL  =  [x(q(L)),x{q(L)-l),---,x(q(L)-(Q-l))} 

=  [x(M),x(M  -  1),  ■  •  ■  ,x(M  -  (Q  -  1))] 

(40) 

There  are  three  cases  to  be  considered:  (1)  g(l)  >  Q,  (2)  1  <  q(  1)  <  Q,  and  (3) 
g(l)  <  1.  It  is  assumed  that  1  <  g(l)  <  Q,  i.e.,  1  <  M/L  <  Q.  The  other  two  cases 
are  different  but  very  similar.  The  input  samples 

x(M),x(M  -  1),  x(M  -  2),  •  •  •  ,x(2),x(l),  x(0),  x(-l),  •  •  • ,  x(g(l)  -  (Q  -  1))  are 
involved  in  the  computation  of  the  first  L  output  samples  ?/(l),  y(2),  ,  y(L).  If  it  is 

assumed  that  the  input  samples  x(n)  start  at  the  sample  x(l),  then  the  samples 
x(0),  ■  •  • ,  x(q(  1)  —  (Q  -  1))  in  the  first  state  variable  buffer  Si  are  not  available  and 
have  to  assume  certain  given  values  (usually  zeros).  Thus  to  compute  y(  1),  the  state 
variable  buffer  S  is  first  initialized  to  be  an  all-zero  vector  and  then  q{  1)  samples 
x(q(l),  x(q(l)  —  1),  x(l)  are  shifted  into  it.  The  samples  in  the  buffer  S  are  then 
weighted  with  the  coefficients  of  the  polyphase  filter  gi(n).  The  most  recent  sample  in 
the  state  variable  buffer  is  x(q(  1)).  To  compute  y( 2),  q( 2)  -  g(l)  new  samples  from 
x(n)  are  shifted  into  S  and  q( 2)  -  q(  1)  old  samples  are  shifted  out  of  S.  The  samples 
in  the  updated  buffer  are  weighted  with  the  coefficients  of  the  polyphase  filter  p2(rc). 
The  most  recent  sample  in  the  state  variable  buffer  has  become  x(g(2)).  In  general,  to 
compute  y(m),  3  <  m  <  L  -  1,  q(m)  -  q(m  -  1)  new  samples  are  shifted  into  the 
buffer  S  and  q{m)  —  q{m  —  1)  old  samples  are  shifted  out  of  the  buffer  S.  The  samples 
in  the  buffer  are  weighted  with  the  coefficients  of  the  polyphase  filter  gm(n).  To 
compute  the  last  sample  y(L)  in  the  first  block  of  L  output  samples, 
q(L)  —  q(L  —  1)  =  M  —  q(L  —  1)  samples  are  shifted  into  and  out  of  the  buffer  S  and 
the  samples  in  the  state  variable  buffer  are  weighted  with  the  coefficients  of  the 
polyphase  filter  g/,(n).  The  most  recent  sample  in  the  buffer  has  become  x(M).  Note 
that  the  first  M  input  samples  x(M)yx(M  -  1),  •  •  • ,  x(l)  have  passed  into  the  buffer  S 
during  the  computation  of  the  first  L  output  samples  y(L),  y(L  -  1),  •  •  • ,  y(  1). 

The  pattern  in  which  the  data  samples  x(n)  are  shifted  into  the  state  variable  buffer  S 
during  the  computation  of  the  next  block  of  L  output  samples 
y(L  +  1),  y(L  +  2),  •  •  • ,  y(2L)  will  now  be  examined.  Denoting  the  contents  of  the 
state  variable  buffer  S  corresponding  to  the  samples  y(L  -f  1),  y(L  +  2),  •  •  • ,  y(2L)  by 
Sl+i,  *  ■  S2 respectively  ,  we  have 

SL+1  =  [x(M  +  q(l)),x(M  +  q(l)  —  !),••  •  ,x(M  +  g(l)  —  (Q  —  1))] 
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Sl+2  = 

Si+3  = 

[x(M  +  q(2)),  x(M  4-  q{2)  —  1),  • 
[x(M  +  <?(3)),  x(M  +  q(3)  -  1),  • 

■  •  ,x(M  +  q(2)  —  (Q  —  1))] 

•  ■  ,x(M  +  9(3)  —  (Q  —  1))] 

S2L  = 

\x{M  +  q(L)),x(M  +  q(L)-  1), 
\x(2M),x(2M  -!),•••,  x(2M  - 

•  •  • ,  x(M  +  q{L)  —  (Q  -  1))1 
(<?-!))] 

(41) 

After  the  computation  of  y(L)  and  before  the  computation  of  y(L  +  1),  the  state 
variable  buffer  S  is 

Sl  =  [x(q(L)),x(q(L)  -  1),  •  •  ■  ML)  -  (Q  -  1))] 

(42)  =  [x(M),x(M  —  1),  •  ’  *  ,x(M  —  {Q  -  1))] 

To  compute  y(L  +  1),  q(l)  new  samples  are  shifted  into  S  and  g(l)  old  samples  are 
shifted  out  of  S.  The  samples  in  the  buffer  are  weighted  with  the  same  set  of 
coefficients  as  in  the  computation  of  y(l) ,  namely,  the  polyphase  filter  gi(n).  To 
compute  y(L  +  2),  q( 2)  —  g(l)  new  samples  are  shifted  into  S  and  out  of  S.  The 
samples  in  the  buffer  are  weighted  with  the  same  set  of  coefficients  as  in  the 
computation  of  y( 2),  namely,  the  polyphase  filter  g2(n).  It  is  now  apparent  that  just  as 
the  polyphase  filters  enter  into  the  computation  of  the  output  samples  y(m) 
periodically,  the  input  data  samples  x(n)  are  shifted  into  the  state  variable  buffer  in  a 
manner  that  also  repeats  itself  periodically.  Thus  the  computation  of  the  output  samples 
y(m)  can  be  done  most  conveniently  on  a  block  by  block  basis.  For  each  block  of  input 
samples  of  length  M,  a  block  of  output  samples  of  length  L  can  be  computed.  This  is 
illustrated  in  the  diagram  of  Fig.  6. 

3.  IMPLEMENTATION  OF  THE  PERIODICALLY 

TIME-VARYING  FIR  FILTER  STRUCTURE  IN  MEX 


Implementation  of  the  periodically  time  varying  FIR  filter  structure  (Eq.  33)  in  some 
software  packages,  including  the  Intel  Signal  Processing  Library  [2],  requires  that  the 
length  of  each  input  data  block  from  the  sequence  x(n)  be  an  integer  multiple  of  the 
downsampling  factor  M.  On  the  one  hand,  this  restriction  considerably  simplifies  the 
indexing  in  the  code  and  reduces  the  complexity  of  implementation.  On  the  other  hand, 
it  also,  to  a  certain  degree,  unnecessarily  limits  the  flexibility  and  usefulness  of  the 
resulting  program,  since  in  certain  applications  the  length  of  the  input  data  blocks  may 
not  be  an  integer  multiple  of  the  downsampling  factor  M.  For  such  applications,  there 
is  a  need  for  an  implementation  of  the  resampling  algorithm  (Eq.  33)  that  imposes  no 
restriction  on  the  length  of  the  input  data  blocks.  The  MATLAB  built-in  function 
UPFIRDN.dll  (or  RESAMPLE.m)  is  designed  for  resampling  only  one  block  of  input 
data.  It  does  not  provide  information  on  the  state  variables  and  cannot  be  used  to 
resample  a  huge  data  file.  To  obtain  a  more  flexible  resampling  program  than  the  ones 
currently  available,  we  have  implemented  the  periodically  time  varying  FIR  filter 
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structure  of  Eq.  33  in  such  a  way  that  it  keeps  track  of  the  state  variables  and  also 
makes  provisions  for  input  blocks  of  arbitrary  length.  The  implementation  is  carried  out 
in  MEX  in  the  MATLAB  environment  with  the  core  of  the  program  written  in  C.  In 
addition  to  utilizing  the  very  efficient  periodically  time- varying  FIR  filter  structure,  we 
also  strive  to  minimize  data  movement  in  our  implementation.  It  turns  out  that  our 
implementation  is  not  only  very  flexible  but  also  works  almost  twice  as  fast  as  the 
Matlab  built-in  function  UPFIRDN.dll  if  the  resampling  ratio  is  non-trivial  (that  is, 

L  >  1  and  M  >  1). 

How  the  program  works  will  now  be  explained.  It  is  assumed  that  the  input  signal 
sequence  is  fed  into  the  computer  memory  on  a  block  by  block  basis.  To  make  the 
program  general  enough,  the  individual  input  data  blocks  are  not  assumed  to  have  equal 
length.  Hence  each  individual  input  data  block  can  be  of  any  length.  Once  a  data  block 
is  fed  into  the  computer  memory,  it  is  resampled  and  the  resampled  data  is  then  stored 
in  the  computer  memory  for  further  processing  or  output  to  a  disk  for  storage.  Since  the 
algorithm  in  Eq.  33  is  much  easier  to  implement  for  an  input  data  block  whose  size  is 
an  integer  multiple  of  the  downsampling  factor  M,  we  divide  each  input  data  block  into 
two  segments.  The  length  of  the  first  segment  is  an  integer  multiple  of  M(in  some 
extreme  cases  it  may  be  empty)  and  the  second  segment  consists  of  less  than  M 
samples  (it  may  also  be  empty).  The  samples  in  the  first  segment  are  processed  in 
blocks  of  M  samples  each  according  to  the  scheme  discussed  in  the  preceding  section 
(see  Fig.  6).  The  remaining  samples  at  the  end  of  the  input  data  block,  which  number 
less  than  M,  are  not  processed.  They  are  returned  in  an  array  called  REMAINDER. 

The  program  will  add  the  samples  in  the  REMAINDER  array  to  the  beginning  of  the 
next  input  data  block  to  form  a  new  elongated  array.  This  new  input  data  block  is  again 
divided  into  two  segments.  The  length  of  the  first  segment  is  an  integer  multiple  of  M 
and  the  second  segment  consists  of  less  than  M  samples.  Again  the  first  segment  is 
processed  and  the  second  segment  is  returned  in  the  REMAINDER  array. 

More  specifically,  the  processing  in  the  program  is  carried  out  in  three  steps  which  are 
depicted  in  the  diagram  in  Fig.  7.  The  first  step  in  the  program  is  to  append  the  old 
remainder  array  of  length  L0  (if  it  is  not  empty)  to  the  beginning  of  the  input  data  array. 
This  step  results  in  an  elongated  input  data  array.  It  is  emphasized  that  this  is 
accomplished  via  pointer  manipulation  and  no  data  are  really  moved  in  the  process. 
After  this  is  done,  the  elongated  input  data  array  is  partitioned  into  two  segments.  The 
length  of  the  first  segment  is  an  integer  multiple  of  M  and  there  are  t  such  blocks  each 
of  size  M  (t  may  be  zero).  The  second  segment  consists  of  less  than  M  samples  (. L\ 
samples).  The  first  segment  is  processed  on  a  block  by  block  basis  as  explained  earlier 
and  the  results  are  returned,  which  make  up  an  output  array  of  length  tL  consisting  of  t 
blocks  each  of  size  L.  The  second  segment  is  returned  as  the  new  remainder  array.  The 
contents  of  the  state  variable  buffer  right  after  processing  the  t  input  blocks  of  size  M 
are  returned  as  the  new  state  variable  buffer  samples. 

In  the  program,  there  is  an  interface  component  (named  the  MexFunction)  that  connects 
the  MATLAB  environment  with  the  actual  data  processing  performed  in  C.  A  large 
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portion  of  indexing  management  is  done  in  this  interface  component.  As  pointed  out  in 
the  preceding  section,  for  each  of  the  three  separate  cases,  q(  1)  >  Q,  or  1  <  g(l)  <  Q 
or  q(l)  <  1,  the  manner  in  which  data  are  shifted  into  and  out  of  the  state  variable 
buffer  is  different  and  the  subscripting  of  the  data  arrays  has  to  be  managed  in  the 
program  accordingly.  To  simplify  the  indexing  and  reduce  the  complexity  of  the  code, 
three  C  functions,  RESAMP1,  RESAMP2  and  RESAMP3,  which,  respectively, 
correspond  to  the  three  separate  cases,  g(l)  >  Q,  or  1  <  q{m)  <  Q ,  or  q(m)  <  1,  are 
written  in  the  program.  Since  the  data  processing  procedure  is  exactly  the  same  for  both 
the  real  and  imaginary  parts  of  the  data  sample  blocks,  the  input  data  sample  blocks  are 
separated  into  the  real  and  imaginary  parts  and  the  C  functions  RESAMP1,  RESAMP2 
and  RESAMP3  are  written  to  accept  double  precision  real  data  arrays  and  to  return 
double  precision  real  data  arrays.  The  C  function  RES AMP2  is  discussed  in  some 
detail  here.  The  other  two  functions  are  very  similar.  The  prototype  for  the  function 
RESAMP2  is  shown  below,  and  the  parameters  for  this  function  are  defined  in  Table  3. 


RESAMP2(  double  Input[  ],  long  InputLen, 
double  Output[  ],  long  OutputLen, 
double  State[  ],  int  StateLen, 
double  Filter[  ],  int  FilterLen, 
double  Remainder^  ],  int  RemainderLen, 
double  Newremainder[  ],  int  NewremainderLen, 
int  L,  int  M  ); 


In  the  C  function  RESAMP2,  the  prototype  filter  array  Filter[  ]  is  rearranged  to  result  in 
a  new  filter  array  named  firfilter[  ],  which  corresponds  to  the  array  of  polyphase  filters 
in  the  upper  row  in  Fig.  6.  Two  arrays,  named  yindxmodLf  ]  and  newyindxmodL[  ] 
respectively,  are  defined  to  hold  the  control  sequences  q(m)  and  p(m),  1  <  m  <  L. 
The  program  completes  the  processing  of  data  blocks  of  size  M  in  the  for  loops, 
returns  the  remainder  samples  and  updates  the  state  variable  samples.  The  complete 
MEX  file  is  named  UPDNC6.C.  After  compilation,  UPDNC6.dll  can  be  invoked 
directly  in  the  MATLAB  environment.  The  usage  of  UPDNC6.dll  and  two  related 
MATLAB  functions  is  discussed  in  the  next  section.  For  details  of  the  C  functions, 

RES  AMP  1,  RESAMP2  and  RESAMP3,  the  reader  is  referred  to  the  file  UPDNC6.C. 

4.  USAGE  OF  THE  FUNCTIONS  UPDNM6.m  AND 
RESAM6.m 


Although  the  program  UPDNC6.dll  can  be  used  directly  in  the  MATLAB  environment, 
we  have  written  another  two  MATLAB  functions  named  UPDNM6.m  and  RESAM6.m 
respectively,  which  are  more  convenient  to  apply  in  certain  situations.  The  function 
UPDNM6.m  calls  UPDNC6.dll  and  the  function  RESAM6.m  in  turn  calls 
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Table  1:  Parameters  for  Resampler  Program 


Parameter 

Description 

Input[  ] 

Input  data  array 

InputLen 

Input  data  array  length,  integer  multiple  of  M 

Output  [  ] 

Output  array 

OutputLen 

Output  array  length,  integer  multiple  of  L 

Statef  ] 

Old  state  variable  buffer  from  earlier  input  data  block 

StateLen 

State  variable  buffer  length  Q  =length(Filter[])/L 

Filter!  ] 

Taps  of  prototype  filter,  h(k),  0  <  k  <  LQ  -  1 

FilterLen 

Length  of  prototype  resampling  filter 

Remainder[  ] 

Old  remainder  array 

RemainderLen 

Length  of  old  remainder 

Newremainderf  ] 

New  remainder  array 

NewremainderLen 

Length  of  new  remainder  array 

L 

Upsampling  factor 

M 

Downsampling  factor 
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UPDNM6.m.  Their  usage  is  explained  here. 


4.1  Usage  of  UPDNM6.m. 

The  MATLAB  function  UPDNM6.m  is  called  with  the  syntax: 

[Y,NEWSTATE,NEWREMAINDER,NEWYINDEX]  = 

(43)  updnm6(X,H,L,M,STATE, REMAINDER, YINDEX) 

where  the  inputs  are: 

1 .  X:  input  signal  with  arbitrary  length,  where  X  can  be  real  or  complex; 

2.  H:  prototype  resampling  filter  designed  by  the  user; 

3.  L:  upsampling  factor,  where  L  must  be  a  positive  integer; 

4.  M:  downsampling  factor,  M  a  positive  integer;  L  and  M  relatively  prime; 

5.  STATE:  old  state  variables,  length(STATE)=length(H)/L; 

6.  REMAINDER:  samples  left  over  from  the  preceding  block  of  input  data; 

7.  YINDEX:  length  of  the  resampled  data  obtained  by  resampling  the  input  data 
before  the  current  input  block  X; 


The  outputs  are: 

1 .  Y:  the  resampled  data  from  processing  all  the  data  blocks  of  size  M  contained  in  X; 

2.  NEWSTATE:  updated  state  variables; 

3.  NEWREMAINDER:  samples  left  over  from  the  current  block  X; 

4.  NEWYINDEX:  length  of  the  resampled  data  obtained  by  resampling  input  data  up 
to  and  including  the  current  block  X; 


When  using  UPDNM6.m,  it  is  necessary  for  the  user  to  supply  the  prototype 
resampling  filter  H.  It  should  be  emphasized  that  there  is  one  restriction  on  the  length 
of  H,  namely,  the  length  of  H  must  be  an  integer  multiple  of  the  upsampling  factor  L.  If 
the  user  does  not  want  to  design  the  prototype  resampling  filter  H,  we  have  written 
another  MATLAB  function  which  designs  a  prototype  resampling  filter  for  the  user. 
This  program  is  named  RESAM6.m.  Its  usage  is  explained  next. 
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4.2  Usage  of  RESAM6.m. 

The  MATLAB  function  UPDNM6.m  is  called  with  the  syntax 

[Y,NEWSTATE,NEWREMAINDER,NEWYINDEX,H] 

(44)  =  resam6(X,L,M,  STATE,  REMAINDER, YINDEX) 

where  the  inputs  are: 

1.  X:  Current  input  signal  with  arbitrary  length; 

2.  L:  Upsapling  rate,  L  must  be  a  positive  integer; 

3.  M:  Downsampling  rate,  M  must  be  a  positive  integer;  L  and  M  must  be  relatively 
prime; 

4.  STATE:  Old  state  variables  after  processing  the  input  data  block  which  precedes 
the  current  input  block  X; 

5.  REMAINDER:  Input  samples  left  over  from  the  input  data  block  which  precedes 
the  current  input  block  X; 


The  outputs  are: 

1 .  Y:  The  resampled  data  from  processing  all  the  input  data  blocks  of  size  M 
contained  in  X; 

2.  NEWSTATE:  updated  state  variables; 

3.  NEWREMAINDER:  input  samples  left  over  from  the  input  block  X; 

4.  NEWYINDEX:  length  of  the  resampled  data  Y. 

5.  H:  FIR  filter  designed  and  used  in  the  resampling  operation; 


The  function  RESAM6.m  does  the  following: 


1.  Creates  a  prototype  resampling  filter  H,  where  H  is  designed  using  the  following 
MATLAB  parameters: 


rp 

r8 

Fa 

f 

a 

dev 

[n,  f0,  a0,  w] 
n 

H 

H 


1  (dB);  — passband  ripple 
60  (dB);  — stopband  ripple 

2  (Hz);  — normalized  sampling  frequency 

[  l/max(L,  M),  1.5/max(L,  M)  ];  — normalized  cutoff  frequencies 
amplitudes 


remezord(/,  a,  dev ,  Fs)\  — filter  order  estimation 
ceil(n/L)  *  L  -  1;  — choosing  filter  length 

L  *  remez(n,  /0,  ao,  w);  — design  filter  with  the  REMEZ  function 
H(:)l 


)];  — deviations  from  the  ideal  response 


[1,  0];  — desired 
l0,_ 


2.  Resample  the  data  block  X  using  the  filter  H  by  calling  the  function  UPDNM6.m. 
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5.  PERFORMANCE  COMPARISON 


To  obtain  an  idea  of  how  our  function  UPDNC6.dll  compares  with  the  MATLAB 
built-in  function  UPFIRDN.dll  in  performance,  we  tested  them  for  various  resampling 
ratios  and  various  input  block  lengths.  The  personal  computer  used  in  the 
simulations  has  a  clock  rate  of  200  MHz.  The  computer  processing  time  (in  seconds)  is 
plotted  as  a  function  of  the  length  of  the  input  block  in  Figs.  8  to  30.  It  can  be  observed 
from  Fig.  8  to  13  that  if  the  downsampling  rate  M  =  1,  our  function  UPDNC6.dll 
outperforms  the  MATLAB  implementation  UPFIRDN.dll  for  shorter  filter  lengths  and 
has  roughly  the  same  performance  as  UPDFIRDN.dll  for  longer  filter  lengths.  If  the 
upsampling  rate  L  =  1,  the  MATLAB  implementation  UPFIRDN.dll  outperforms  our 
implementation  UPDNC6.dll  for  shorter  filter  lengths  but  has  roughly  the  same 
performance  as  ours  for  longer  filter  lengths.  This  is  demonstrated  in  Figs.  14  to  19.  If 
the  resampling  ratio  is  non-trivial,  that  is,  if  L  >  1  and  M  >  1,  then  our 
implementation  UPDNC6.dll  outperforms  the  MATLAB  implementation  UPFIRDN.dll 
by  a  factor  of  about  2  for  all  filter  lengths.  This  is  demonstrated  in  Figs.  20  to  30. 

6.  CONCLUDING  REMARKS 


We  have  successfully  implemented  in  software  the  periodically  time- varying  FIR  filter 
structure  in  digital  sampling  rate  conversion  and  obtained  a  very  flexible  resampling 
function  UPDNC6.dll.  This  program  outperforms  the  MATLAB  built-in  function 
UPFIRDN.dll  by  a  factor  of  about  2  for  non-trivial  resampling  ratios  (that  is,  L  and 
M  are  relatively  prime  and  L  >  1  and  M  >  1). 
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Figure  8:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  5,  downsampling 
factor  M  -  1,  and  filter  length  =  20. 
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Figure  9:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  5,  downsampling 
factor  M  —  1,  and  filter  length  =  40. 
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Figure  10:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  25,  downsampling 
factor  M  =  1,  and  filter  length  -  160. 
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Figure  11:  Comparison  of  MAT  LAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  25,  downsampling 
factor  M  =  1,  and  filter  length  ~  125. 
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Figure  12:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  25,  downsampling 
factor  M  —  1 ,  and  filter  length  =  250. 
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Figure  13:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  ~  25,  downsampling 
factor  M  =  1,  and  filter  length  =  1000. 
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Figure  14:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  1,  downsampling 
factor  M  —  5,  and  filter  length  =  20. 
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Figure  15:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  1,  downsampling 
factor  M  —  5,  and  filter  length  =  1 00. 
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Figure  16:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  1,  downsampling 
factor  M  =  5,  and  filter  length  =  500. 
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Figure  1 7:  Comparison  of  MATLAB  function  UPFiRDN  and  UPDNC6,  for  upsampling  factor  L  -  1 ,  downsampling 
factor  M  —  25,  and  filter  length  =  200. 
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Figure  18:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  1,  downsampling 
factor  M  —  25,  and  filter  length  =  600. 
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Number  of  Samples  x1Qs 


Figure  20:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  —  5,  downsampling 
factor  M  =  4,  and  filter  length  =  50. 
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Figure  21:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  5,  downsampling 
factor  M  =  4,  and  filter  length  =  100. 
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Figure  22:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  —  5,  downsampling 
factor  M  =  4,  and  filter  length  =  200. 
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Figure  23:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  25,  downsampling 
factor  M  =  24,  and  filter  length  =  125. 
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Figure  26:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  25,  downsampling 
factor  M  =  24,  and  filter  length  =  3000. 
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Figure  27:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  24,  downsampling 
factor  M  =  25,  and  filter  length  =  192. 
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Figure  28:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampiing  factor  L  =  24,  downsampling 
factor  M  =  25,  and  filter  length  =  480. 
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Figure  29:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  24,  downsampling 
factor  M  =  25,  and  filter  length  =  960. 
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Figure  30:  Comparison  of  MATLAB  function  UPFIRDN  and  UPDNC6,  for  upsampling  factor  L  =  24,  downsampling 
factor  M  —  25,  and  filter  length  ~  2400. 
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